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We derive the Wick theorem for the q-Exponential distribution. We use the theorem 
to derive a numerical algorithm for finding parameters of the correlation matrix of q- 
Exponentialy distributed random variables given empirical spectral moments of the time 
series. 
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1. Introduction 



q- Exponential distributions (q- Exponentials defined in the next section) are used[l] 
in modeling distributions of stocks. The q- Exponentials possess two desirable fca- 
' tures namely exhibit power law tails in the high end of the distribution and tend 

toward a Gaussian when g = 1 + l/i"T — > 1. These facts provide enough motiva- 
tion to derive a variant of a Wick theorem for linear combinations of independent 
identically distributed (iid) q- Exponentials (correlated q- Exponentials) . 

Recent growth of interest in applications of physics to economics and to the 
theory of finance (econophysics) yields the Wick theorem useful for the following 
purposes: 



(1) Deriving exact relations between spectra of eigenvalues of financial covariance 
matrices related to g-exponentially distributed time series and of estimators of 
covariance matrices. This will be a generalization of existing exact relations [8] 
for Gaussian distributed time series. 

(2) Deriving algorithms [6] for optimizing a portfolio (minimizing the variance 
and/or the higher moments of a portfolio subject to a presumed return from a 
portfolio) of stocks whose time series are g-exponentially distributed. 

(3) Verifying empirical findings [5] regarding a very fast time decay of the autocorre- 
lation function in time series and regarding a power-law time dependence of the 
autocorrelation function of the volatility; investigating many point correlations 
in financial time series. 
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2. The Wick theorem 

The objective of this section is to derive the Fourier transform of N correlated 
g-Exponentially distributed variables. 

A random variable X is q-Exponcntially distributed, (X <~ q — Exp) if q = 
1 + 1/K and the probability density function (pdf) D K (x) of X reads: 



D K {x) 



N k 



(2.1) 



V^l K 

wh ere e z K := ( 1 - z/K)- R , N k = r(K)/(VK ■ T(K - 1/2)) and a K = 
a^(K - 3/2)/K. 

The pdf (2.1) has (equation (20.52) in [1]) is a continuous superposition of 
Gaussians and has a following integral representation: 



(2.2) 



where D := ^/2tt/(D + 1) and a D = <r^(D - 2)/(D + 1) = u K and the integral 
runs over the D-dimensional space • • • , £d) with: 

D = 2K - 1 (2.3) 

Wc define correlated g-Exponentials as linear combinations of iid g-Exponentials. 
Therefore the joint pdf p x (x) of NT correlated (/-Exponentials 

N T 

^ = EE°tfe ( 2 - 4 ) 

j=i t=i 

where i = 1, . . . , N, t = 1, . . . , T and Yj_g <~ g — Exp are iid, reads: 

- AT T 



= (det£)" 1 [] £ (O" 1 )^,*) 



(2.5) 
(2.6) 



»=1,...,JV 
t=l,...,T 



j=l,...,JV 

e=i,...,r 



Here we call O := |o| ; f } j _ ;[ AT a ro ^ at ^ on t cnsor - We take k := {fcj,t} and 

i = l,...,T 

we calculate the Fourier transform n{k) = T x [p x ] (k) := J r jvt d NT xp x (x)e lkx . It 
reads: 



(2) DArT detO) K (fc) = (V DNT dctO) 



O) I d* 



xp x {x)e lkS = 



( 



n 



»=1,...,JV 
t=l T 



n 



-« l ? t (£'+l)/2 



i=l,... 
>t=l,...,T 



rri V 



jNT - 

a xexp 



2^ 



. j=l,...,N 
\9=1,...,T 
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/ -/ 



n d ^ 

i=l,...,N 
t=l,...,T 



i=l,..., N 
\t=l,...,T 



d NT x exp \-^x T ■ C-\i) ■ x + ikx\ 

R NT { 2 — J 



where we introduced (7(£) := O - D O 7 (a correlation tensor) that is related to the 
rotation tensor O and to a diagonal tensor Df^ = SijSt.ecro/^ff The transposition 
T operation is defined as (0 T ){f := Oj'* e . This means that: 



p=l,...,JV 
A=1,...,T 



(2.9) 



The last integral on the right hand side in (2.8) is evaluated by "completing to 
a square" and it reads: 



(2.)^/ 2 (detfC- 1 ^)- 1 ' 2 ) .expl-ip-g^} 
= (2.)^ (detfi) ^ T f ft ft tf) • cxp 



Therefore the Fourier transform (2.7) reads: 



(2.10) 
(2.11) 



K(jfe) := 

(27T) JVT / 2 (7^ T 



qDNT 

1 



n 



e 



d e (D+l)/2 



j=l,...,N 
9=1, ...,T 



rr\DNT I 



/ n ^(e-^ e(Z?+1)/2 )exp{-ifc r 0a4 



(2.12) 

1 1 exp{-lfc T <2(04 



j=l,...,N 
9=1, ...,T 



■{-ip-iaat}) 



= <«p(-2*'c«)*! 

The weight w(£) used in the average ()|* m (2-14) reads: 



j=l,...,N 
9=1,. ..,T 



(2.13) 
(2.14) 



(2.15) 



where A/" := ((£> + l) D ' 2 )/{2 D / 2 - l T{D /2)). 

In (2.15) we expressed the integral over £j g in radial coordinates J RD eP£ = 

2W^/ 2 /(rp/2))/ 00 ^- 1 ^. 

Now we differentiate (2-14) 2k times with re- 
spect to variables fcj(i),9(i), • • • , fcj(2fe),9(2fc) an d evaluate the result at fc = 0. We 
get: 
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' • ■ • ' x j(2k),e(2k) 



) = 



d 



i 2k dk j{1)M i) 



dk 



■j(2k)M(2k) 



( 



X 



(-l) fe dk o[1)fi (i) " dk ]{2 k),e(2k) \2 k k\ 



i,j=l,...,N 
i,6=l,...,T 



E 



o...oc 



(fe) \g=l 



A^(c<«>),0(4«>),a\ 

llS-(cW),fl(cW)«7 
9=1 / pc 



E 

<t=c< 1 )o...oc(' c ' \<?=i 



n 



lfc=0 



(2.16) 



(2.17) 



(2.18) 



1 1 1=1, ...,N I 




o 



,p(q)Mq) 



o 



m=l,...,T 

n ( (£±i) r « p g-;:<'-'">,t i>w 

m=l,...,T 



r(D/2) 



where the sum in (2.18) runs over all 2fc-permutations a that are composed exclu- 
sively of cycles := (c^\ c^) for i = 1, . . . , k of lengths two. Since the functions 
Cf 'f are symmetric with respect to exchanging the above and lower pairs and mul- 
tiplication is commutative the whole of (2k)\ terms in the sum in (2.17) decomposes 
into (2fc)!/(2 fc fc!) = (2k— 1)!! distinct terms who occur (2 k k\) times each. In this way 
the factor 2 k k\ in the denominator in (2.17) cancels out. In (2.19) we used the defi- 
nition (2.9) of the correlation tensor and we introduced new indices p(q) =1,...,N 
and X(q) = 1, . . . , T Since the average over £ in (2.18) consists in performing inte- 
grals of the kind Af J Rd Yl q =i Oq) x^^iO^ ^ is readily seen that the result (2.19) 
is expressed via a number v(l,m) that depends on the sequence {p(q),X(q)} and 
that is equal to the multiplicity of the pair (l,m) in that sequence: 



m) = # {pairs (p q , X q ) such that p q = I and \ q — m} 



(2.21) 



3. Covariances their estimators and noise filtering 

In this section we discuss definitions of averages over ensembles of stochastic vari- 
ables (resolvents) whose properties may be compared to measured properties of 
financial time series. 

Definition 3.1. The resolvent function is a complex function G(z) such that: 

Im[G(z + ie)]=S e (z), e>0 (3.1) 
where 6 e (z) is a representation of the delta function. 
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The function G(z) is used for finding the density of eigenvalues Pa(A) of the 
covariance (cQQj.j) where 



1 T 



(3.2) 



t=l 

and the average is over the random ensemble X i>t . We have: 

N N 



p A (X) = lim e ^o J2 <^ A - A *) = lim ^° J2 Im[G(A - A 4 + ie)] (3.3) 
»=i i=i 
= lim e ^oTr [im [G(A1 -g + ie)]] (3.4) 



Lemma 3.1. T/ie resolvent function G(z) has res z= iG(z) — l/2n. Here res de- 
notes a residue. 



Proof. The function G(z) can be expanded in a Laurent series around z = 0. This 
means that 



3 fe > G(z) = Y / ^+g(z) 



3 = 1 



(3.5) 



where the function Q(z) is analytic. 



We need to prove that Im 
the non-analytic term in (3.5). 

Im 



/ G(x + l e) 



1 for e > 0. We consider at first 



1 



(x + ie)i 



dx 



1 3 = 1 
j > 1 



(3.6) 



where the above equality in (3.6) is straightforward and the lower equality is derived 
by means of the Cauchy integral theorem. Now we analyse the analytic term. For 
R > we compute 



R+ie 



j g{z + ie)dz= J g(z)dz= j G(z)dz + j \g{R + ie) -Q{-R + ie)]dt, (3.7) 

-R -R+it -R 

where the last equality in (3.7) follows from the application of the Cauchy integral 
theorem to a contour consisting of four intervals [— R,R], [R,R+ie], [R + ie, — R+ie] 
and [— R + le, —R\. From the last expression on the right hand side in (3.7) we see 
that the imaginary part of that integral is 0(e) and hence disappears when e — > 0. 
This fact and (3.6) suffices to finish the proof. □ 



The class of functions G(z) becomes narrowed down subject to the following 
condition: 
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Lemma 3.2. The resolvent function G(z) = 1/z, if and only if 

V neN — ^(-1)" = (G(z))" +1 (3.8) 

The necessity follows in a straightforward manner from substituting G{z) = 1/z 
into (3.8). To prove the sufficiency we take x € R, multiply both sides of (3.8) by 
x n and sum over n = 0, 1, . . . , 00. Since the left hand side is the Taylor expansion 
of G(z — x) around z and the right hand side form a geometric series we obtain a 
functional equation: 

G(z-x) = - gg- = ^J— (3 .9) 



1 - xG(z) 



X 



We substitute x = (1 — a)z into (3.9) for some a e R+ and get: 
1 1 



G{az) G{z) 



+ (a-l)z (3.10) 

1 

W) 

= (a-Vjz—^- = az (3.12) 

a 

where in (3.11) we have iterated the equation (3.10) and we used the fact that 
1/G(0) = 0. Finally in (3.12) we summed a geometric series and completed the 
proof. 



4. Computation of the expansion of the resolvent 

We calculate a function 

g(*) = (G(Vl-^-^ T )) (4.1) 

where the average is performed over random variables (2Qi,e ~ q — Exp that are 
correlated in i and in time 9 and G(z) is a resolvent function. The function g(z) 

is termed the resolvent. For the purpose of the calculation we fix ^ € M. D , we 
project the joint pdf p^(x) onto Gaussians Normal(0,<7fj/£? t ) using the integral 
representation (2.2), average over the Gaussians and obtain 

(X ijB X$ J ) = Ctf'®= £ ^#-Ojf (4.2) 

p=l,...,JV ^p. a 
A=1,...,T 

Note that the two-point correlation function (4.2) does not factorize into functions 
depending on times 9,9' and on i, j only. It factorises for |& )t | = 1 if the rotation 
tensor 0^' t A factorises, ie Of' t A = /fT t A . For generic value of £, however, correlations 
in i and in time are coupled with each other. In the following we sum over repeated 
indices (use Einstein's summation convention). The N-type indices and the T-type 
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indices (running from one to N and to T) are denoted by Latin and by Greek letters 
respectively. We expand the resolvent function G() in a series around z. We get: 



n=0 



= Y1 T n n \ ((^•^ T )r(i), 4 ( 2n+ i)/^(2«+l)j 
n=0 

^ r"n!^ _1 ^ 5i ' i(1) (A ( X '( 2 J'- 1 ).€( 2 J'- 1 ) X l(2j-i),i(2j+i; 



) I / <5»(2n+l),j 



(4.3) 

(4.4) 
(4.5) 



n=0 



\j=l 



n=0 



\j=l 



n=0 (7 \ \ j'=l / j'=l ^ ' j'=l ^ ' / 



In (4.4) we introduced two new N-type indices and i(2n + 1); in (4.5) 

we expanded the n-th power (X ■ X T ) .... , and we introduced indices 

1 ^ = / «(l),i(2n+l) 

i(l),i(3), . . . , i(2n — 1) together with £(1), . . . , £(2n — 1). In (4.6) we inserted n 
T-type Kronecker delta functions between the ordered X, X T pairs and n — 1 N- 
type delta functions between the pairs X T , X. This resulted in introducing indices 
£(2), . . . , £(2n) and indices i(2), . . . , i(2n). 

Finally in (4.7) we made use of the Wick theorem (2.18) (the sum over a runs 
over all 2n-permutations composed entirely of cycles of length two), of Lemmas 
3.1 and 3.2 and we eliminated the index i(2n +1). Note that the sum over a 
contains (2n — 1)!! terms each containing An variable indices i(l), . . . ,i(2n) and 
£(1), . . . , £(2n) and two fixed indices i and j. 

Now we construct a pictorial representation (Feynmann diagrams) [2,3,4,7] of 
terms in the sum (4.7) according to a following recipe. The i(. . .) indices are denoted 
by bullets •, and the £(. . .) indices are denoted by open circles o. The factors in the 
product on the right hand side in (4.7) are assigned to graphs as follows: 
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Factor Graph 




I*9(l),0(2) 



lV d i(l),i(2) 

In this way the sum (4.7) can be represented as a sum over graphs that are con- 
structed from building blocks (4) in such a way that for every node, except for the 
nodes i and j (external nodes), there are exactly two edges abutting at it (this 
follows from the Einstein's summation convention). All graphs contributing to the 
second and to the third order of the expansion are listed, together with their weights, 
in Figs. 3 and 4. 

Each graph consists of a number of closed solid, a number of closed dashed loops 
and of a solid line that starts at index i and ends at j. Since a closed loop corresponds 
to a contraction (setting two indices of the tensor equal and summing over them) 
°f a tensor that is constructed by multiplying the weights of the 

edges of the graph, the weight of the closed loop is proportional to N and to T for 
solid and a dashed closed loops respectively. In other words the weight of a closed 
solid (dashed) loop is equal to a trace of a N x N or (T x T) matrix and thus 
proportional to N or (T). In the following we assume that N/T = r is fixed and 
investigate the expansion (4.7) in the limit N — > oo. Only such graphs contribute to 
the expansion whose number of loops (cither dashed or solid) is equal to the order 
of the expansion n (planar graphs). Graphs that consist of intersecting lines, like 
the second graph from the top in Fig. 3 are negligible in the limit N — > oo since 
their weight is inverse proportional to a certain power of N. 

We define one-line irreducible graphs as graphs that cannot be split into two 
distinct graphs by cutting a certain edge. The usefulness of this definition follows 
from the fact that the weight of a graph (the term in the sum (4.7) for fixed n and 
fixed a) factorizes into a product of weights corresponding to one-line irreducible 
components. Note that this would not be the case if the weight depended explicitly 
on n (e.g. the function G(z) satisfied Lemma 3.1 but not Lemma 3.2) or if it 
depended on some function of indices i{. . .) or £(...)). Denoting by S = j^i}. 

the sum of weights of all one-line irreducible graphs (self-energy) with external 
nodes i and j we realize from (4.7) and from Fig.l that the resolvent is a sum of a 
geometric series in self-energy: 
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. # + . . + . £Z\ . {Z^ . + ... 

i j i ii i 2 j i ii i 2 i 3 "4 i 6 j 

Fig. 1. Computation of the resolvent as a geometric sum of self-energy graphs. 

£(*) = ^ + ^1 + ^sif + • • • = ( zl ~ D _1 ( 4 - 8 ) 

This result is quite general, ie it holds also if non-planar graphs are taken into 
account. However, the self-energy can only be calculated in the case when non- 
planar graphs are neglected as we see in Fig. 2 and in equation (4.9). 

1 1 6 6 1 6 9 

= Ct£((Tl-§)-Xl (4.9) 

where B 8 ^ := Cp^g(z) Piq . Denoting i = i(0) and j = i(n) we write the result in a 
compact way: 

oo -. /n—1 \ 

n=0 \q=0 / 

oo f<j, e ( m + 1 ) (A m 

4 = E iM tU II c%r^)^)m-^m) (4.1D 

m— q— 1 

and we make following comments: 

(1) The results (4.10) and (4.11) are to be understood as follows. We fix £ = 

• ■ • ! £n,t) € M. NT ' D and for a given rotation tensor O and ou we compute 
the correlation tensor C(£) from (2.9), we insert the self-energy from equation 
(4.11) into (4.10) and we iterate the result until convergence, in the expansion 
in 1/z to a given order, is obtained. Then we average the result over £ e 
H NT ' D according to the rule (2.15). Since the weight in (2.15) factorises and 
the resolvent as a function of £ is a polynomial of inverse powers of £j.t the 
weighting can be done analytically. 

(2) The correlation tensor C_ is not a tensor product of i- and time t- dependent 
matrices (does not factorize) (as noted in the paragraph under equation (4.2)) 
even if the rotation tensor O factorises. 

(3) The result (4.11) is only valid in the limit TV — > oo. For finite values of N there 
will be corrections proportional to inverse powers of N, corrections resulting 
from non-planar graphs (see Figs. 3 and 4). 
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Fig. 2. Calculation of the self energy as a geometric sum of the resolvent contributions. 



5. The direct and inverse problems and the moment expansion of 
the resolvent 

The purpose of this section is to make connection between quantities that are mea- 
sured from financial time series, namely estimators of correlations c := cQQij 
(definition (3.2)), and between the underlying rotation tensor Off and the vari- 
ance (J d of the q-Exponential distribution. We define spectral moments m n of the 
estimator of correlations c as traces of powers of the estimator: 



N 



Tr[c n ] 



(5.1) 



for n = 0, 1, 2 . . .. Having done that we readily see from the definition of the resolvent 
(4.1) that the trace of the resolvent t(z) :— -^-Tr g(z) is a generating function of 

the spectral moments t(z) = X^^Lo m n/z n+1 

Following the terminology from [8] we define two problems the direct and the 
indirect one. The direct problem consists in computing the spectral moments from 
the rotation tensor 0\f and the variance <jd- The indirect problem, which is more 
interesting from the point of view of applications in quantitative finance, is defined 
as determining the correlation tensor and the variance from the sole knowledge of the 
spectral moments. To what extent it is possible, what additional assumptions about 
the structure of the tensor have to be made before deriving a numerical algorithm 
and what is the error estimate in the algorithm will be discussed in future work. 



5.1. Finding the resolvent 



We solve equations (4.10) and (4.11) for g(z)ij according to the recipe 
in point (1) at the end of section 4 and obtain following results: 
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Order n The resolvent g^ n \z)ij 














*4 




(5.2) 


1 






u 

z 




f ?(f c ''*') 




(5.3) 


2 


1 . 

z 


1 

- 


T \ / 


+ 


z 3 7^2 V i' 61 ! 


■ s-lkfil s-lj,$2 


(5.4) 




1 r 1 

z z z 


1 

T 


+ 


1 

" «3 


'■ / /~ijfo/~ihfi2 , 
y2 + 




(5.5) 



O 11 

C k J:cZ\Ci}f 3 + C^CllfC^) (5.6) 
Comparing the coefficients of the expansion (5.6) in powers of l/z with the spec- 
tral moments m„ of the estimator of the correlations we obtain a set of non-linear 
equations that relate certain contractions of the correlation tensor to the spectral 
moments. If we assumed that the correlation tensor factorized, which is not the 
case as we discussed in point (1) in section 4, then we would have obtained equa- 
tions (34) from [8] , equations that relate the spectral moments to moments of the 
underlying correlation matrices both in i and in time. In our case the relations are 
averaged over £ and will be related to some contractions of the rotation tensor O. 
Before proceeding further we note that the relations solve the direct problem but 
they do not provide enough information to solve the indirect problem. 

6. Conclusions 

We have derived a variant of the Wick theorem that expresses the many-point corre- 
lation function of g-Exponentialy distributed random variables through two-point 
correlation functions. This theorem will be used for solving the indirect problem 
in quantitative finance, ie for determining the correlations of time series from the 
knowledge of the spectral moments of the estimator of covariance. 
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